Amplitude dependent frequency, desynchronization, and stabilization in noisy 

metapopulation dynamics 
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The enigmatic stability of population oscillations within ecological systems has remained an open 
theoretical issue for over seven decades. The key models for deterministic description of popula- 
tion dynamics - the Lotka-Volterra prey-predator model and the Nicholson-Bailey host-parasitoid 
system - fail to support an attractive manifold. Accordingly, in any ecological system described by 
these models, species must undergo extinction. Dozens of studies regarding this ecological paradox 
have revealed the important role played by spatial migration and noise, but have yet to pinpoint 
the generic stabilizing process. This underlying mechanism is presented and analyzed here in the 
framework of two interacting species free to migrate between two spatial patches. We show that the 
combined effects of migration and noise cannot account for the stabilization. The missing ingredi- 
ent is the dependence of the oscillations' frequency upon their amplitude; with that, noise-induced 
differences between patches are amplified due to the frequency gradient. Migration among desyn- 
chronized regions then stabilizes a "soft" limit cycle in the vicinity of the homogenous manifold. 
A simple model of diffusively coupled oscillators allows the derivation of quantitative results, like 
the functional dependence of the desynchronization upon diffusion strength and frequency differ- 
ences. Surprisingly, the oscillations' amplitude is shown to be (almost) noise independent. The 
results are compared with a numerical integration of the marginally stable Lotka-Volterra equa- 
tions. An unstable system is extinction-prone for small noise, but stabilizes at larger noise intensity. 
The methodology presented here may be applied to a wide class of spatially-extended ecological 
problems. Indeed, this understanding will likely aid in designing strategies for preventing spatial 
synchronization responsible for species extinction. 



The idea that species populations fluctuate in time has 
been well known since the early days of history. Ancient 
day naturalists, like Herodotus and Cicero, perceived the 
persistence of prey species in the face of adversity as a 
manifestation of divine power and the creator's design 0. 
In modern times, the mathematical description of prey- 
predator interacting populations was given, using deter- 
ministic, continuous time partial differential equations, 
by Lotka and Volterra 0, 0, 01, the analogous model 
with discrete time step was introduced for a parasitoid- 
host system by Nicholson and Bailey Both mod- 
els allow, essentially, for population oscillations around 
a steady state. As pointed out by Nicholson |(J, these 
oscillations are an intrinsic property of interacting popu- 
lations. If the density of the host, say, is above its steady 
value, it will be reduced by the parasite. However, when 
the host reaches its steady density, the density of para- 
sites will be above its steady value. " Consequently, there 
are more than sufficient parasites to destroy the surplus 
hosts, so the host density is still further reduced in the 
following generation ... Clearly, then, the densities of the 
interacting animals should oscillate around their steady 
value" Oscillations in populations and metapopu- 
lations have been observed in many field studies |^ and 
even in controlled experiments The stabilization of 

such oscillations is considered to be a major factor affect- 
ing species conservation and ecological balance 0, ^} . 

Lotka, Volterra and Nicholson recognized that the os- 
cillations described by their models are not stable 0, |(| . 
The Nicholson-Bailey map admits an unstable steady 
state where the amplitude of oscillations grows expo- 



nentially with time; for the Lotka-Volterra system, the 
fixed point is marginally stable, rendering the system 
extinction-prone for any noise amplitude (See Appendix 
1). Indeed, experimental and theoretical studies of both 
systems reveal that the oscillations increase in size until 
one of the species becomes extinct 0, 0, • 

Nicholson Q was perhaps the first propose the idea 
of migration induced stabilization. Although on a sin- 
gle patch, the oscillation amplitude grows in time and 
the system is driven to extinction, desynchronization be- 
tween weakly coupled spatial patches, together with the 
effect of migration, leads to the appearance of spatial pat- 
terns and stabilizes the global populations. This seminal 
idea has been examined in many studies and the main 
results, summarized in a recent review article |14|. are as 
follows: 

• For any network of N patches, if the migration be- 
tween patches is symmetric or almost symmetric 
(i.e., the diffusion of the prey and the predator are, 
more or less, the same), there is no diffusion in- 
duced instability, and the homogenous manifold is 
stable 0,0,03- Thus, the effect of migration 
alone does not cure the instability problem. Diffu- 
sion induced instability may occur if the migration 
rate of the predator is much smaller than that of 
the prey 0] , or in a case where the reaction param- 
eters vary on different spatial patches [TiJ |2(1 |2lJ . 

• Numerical simulations demonstrate, indeed, that 
Lotka-Volterra or even Nicholson Bailey dynamics 
on spatial domains are much more stable [2^, [2^, 
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124. |25j . In general these numerical experiments in- 
volve some sort of noise, like the intrinsic noise due 
to the stochasticity associated with discrete indi- 
viduals (individual-based models), numerical noise 
etc. 

• Accordingly, there is broad agreement that the 
combined effect of noise and diffusion is a necessary 
precondition for population stabilization. However, 
up until now the qualitative nature of the underly- 
ing mechanism has remained obscure, and no theo- 
retical framework that allows for quantitative pre- 
diction has been presented. 

This theoretical gap may be addressed using a toy 
model for coupled oscillators (see Appendix B). The 
main new ingredient emphasized by the proposed model 
is the dependence of frequency on the oscillation ampli- 
tude, reflected by the gradient of the angular velocity 
along the radius oj'(r). The instability induces desyn- 
chronization iff the small, noise-induced, differences be- 
tween patches are amplified by the frequency gradient 
such that the "desynchronization parameter" (0 2 ) ac- 
quires a finite value, leading to "restoring force" toward 
the origin of the homogenous manifold. 

The coupled oscillators model provides the basic the- 
oretical framework with which to explain the emergence 
of an attractive manifold. In the context of the realistic, 
Lotka-Volterra model, the lifetime of the system (time 
until extinction of one of the species) is controlled by 
that manifold. A two-patch system is then described by: 

^ = -^ai + Aidi&i + D a (a 2 - cii) 
do, 

— = -fia 2 + Aia 2 6 2 + D a (ai - a 2 ) (1) 
db\ 

— = ah - A 2 ai6i + D b (b 2 - &i) 
db 2 

— = ab 2 - \ 2 a 2 b 2 + D b (b x - b 2 ). 
ot 

The invariant manifold is the two dimensional subspace 
cii = a 2 , b\ = b 2 . The time evolution of that system, 
with an additive noise, equal diffusivities, D a = D b = D, 
and symmetric reaction rates (See Appendix A) is ob- 
tained through Euler integration. In the limit D = 
the patches are unconnected; thus, starting from the ho- 
mogenous fixed point ax = a 2 = b\ = b 2 = 1, the single 
patch situation (Appendix A) still holds and the system 
hits the absorbing walls after a characteristic, noise de- 
pendent, time. In the opposite limit, D = oo, the system 
sticks to the invariant manifold and acts like a single 
patch (with modified noise and interaction parameters), 
performing again a random walk in the invariant man- 
ifold. However, between these two extremes, there is a 
region where the combined effect of diffusion and noise 
stabilizes a finite region within the invariant manifold. 

The finite nature of the two-patch system ensures that 
it will reach the absorbing state (hit the walls) as t — > oo. 
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FIG. 1: The typical persistence time as a function of the dif- 
fusion rate for different levels of noise. The values of r were 
gathered from survival probability plots (like those in Figure 
3) and are displayed here for the two-patch system. One sees 
that the value of r grows very rapidly (even faster than expo- 
nentially) with the migration rate for small diffusion values, 
and decays with D for large diffusivities. Data is shown for 
different noise intensities A = 0.3 (triangles), 0.5 (squares) 
and 1.0 (circles). 



However, if there is an attractive region in the four di- 
mensional phase space, the death of the system is caused 
by rare events and the typical death times grow consid- 
erably. For any noisy two-patch system Q(t) histograms 
(like those shown for a single patch in Figure 4) were used 
to extract the typical decay time t(D, A) by fitting its 
tail to exponential decay exp(—t / V) . In Figure 1, t(D, A) 
is plotted against D for different noise amplitudes, and 
is shown to increase (faster than exponentially) with D 
(1/D) as D approaches zero (infinity). Evidently, for in- 
termediate diffusivities, an attractive manifold appears 
in phase space, with a Lyapunov exponent that grows 
(faster than linearly) with the diffusion constant. 

On the other hand, at the vicinity of the homogenous 
fixed point, the dynamic is similar to a single patch dy- 
namic. The square of the average distance from the fixed 
point grows linearly with time at the beginning, with a 
slope that depends on the noise amplitude, as expected 
for the random walk in the invariant manifold scenario. 
For "intermediate" migration (e.g., D = 0.01), the aver- 
age distance from the origin saturates, while the chance 
to find the system at large H becomes exponentially 
small, as illustrated in Figure 2. In analogy with the 
results of the toy model, the flow toward the center is 
correlated with the desynchronization, leading to stabi- 
lization of a soft limit cycle at finite H . As predicted, 
while the width of the </> 2 distribution depends strongly 
on the noise amplitude, the oscillation amplitude is al- 
most noise independent. 

In conclusion, we suggest a novel solution to a long- 
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FIG. 2: The average AH at an elementary time step (0.001 of 
a unit time) as a function of the angle (p between the patches. 
While a simple phase space random walk yields on average 
positive AH, this property is shown here to hold only for 
small <f>. At larger angles, the diffusion between patches forces 
the system toward the center and the average AH becomes 
negative. Results are shown for A = 0.1 (full line) and A = 1 
(dashed line). The inset shows the probability distribution 
function for H at these two noise levels. 

standing conundrum: the stabilization of a noisy unsta- 
ble dynamical system on spatial domains. The basic fea- 
ture that leads to stabilization is the dependence of an- 
gular velocity on phase space coordinates. This depen- 
dence allows the noise to desynchronize spatially-coupled 
patches, and then migration decreases concentration gra- 
dients and leads to stabilization of a soft limit cycle close 
to the homogenous manifold. Clearly, owing to the phase 
space confinement of the trajectories, the exact nature 
of the noise becomes quite irrelevant, and thus our con- 
clusions are also applicable to multiplicative noise (indi- 
vidual based dynamics). Furthermore, the effect of rare 
events - the trapping of the dynamics in the absorbing 
state for any finite system - should vanish in the ther- 
modynamic limit, and thus one may expect a parametric 
dependent phase transition in that limit. 
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APPENDIX A: THE NOISY LOTKA-VOLTERRA 
MODEL 

The Lotka-Volterra predator-prey system is a paradig- 
matic model for oscillations in population dynamics 
Hi El LI It describes the time evolution of two interact- 
ing populations: a prey (6) population that grows with 



a constant birth rate a in the absence of a predator (the 
energy resources consumed by the prey are assumed to 
be inexhaustible), while the predator population (a) de- 
cays (with death rate fj,) , in the absence of a prey. Upon 
encounter, the predator may consume the prey with a 
certain probability. Following a consumption event, the 
predator population grows and the prey population de- 
creases. For a well-mixed population, the corresponding 
partial differential equations are: 

NT I A -, S 

— = — fia + Aidb (Al) 

— = ab - \ 2 ab 
ot 

where Ai and A2 are the relative increase (decrease) of 
the predator (prey) populations due to the interaction 
between species, correspondingly. 

The system admits two unstable fixed points: the ab- 
sorbing state a = b = and the state a = 0, b = 
00. There is one marginally stable fixed point at a = 
a/ A2, b = fx/Xi- Local stability analysis yields the 
eigenvalues ±iy/jl<j for the stability matrix. Moreover, 
even beyond the linear regime there is neither conver- 
gence nor repulsion. Using logarithmic variables z = 
ln(a), q = ln(b) eqs. IjAlll take the canonical form 
z = dH/dq, q = —dH/dz, where the conserved quantity 
H (in the ab representation) is: 

H = Ai& + A 2 a - p. ln(a) — a ln(b). (A2) 

The phase space, thus, is segregated into a collection of 
nested one-dimensional trajectories, where each one is 
characterized by a different value of H, as illustrated in 
Figure 3. Given a line connecting the fixed point to one 
of the "walls" (e.g., the dashed line in the phase space 
portrait, Figure 3), H is a monotonic function on that 
line, taking its minimum H m in at the marginally stable 
fixed point (center) and diverging on the wall. Without 
loss of generality, we employ hereon the symmetric pa- 
rameters fj, = <t = Ai = A2 = 1. The corresponding phase 
space, along with the dependence of H on the distance 
from the center and a plot of the oscillation period vs. 
H, are represented in Figure 3). 

Given the integrability of that system, the effect of 
noise is quite trivial: if a and b randomly fluctuate in 
time (e.g., by adding or subtracting small amounts of 
population during each time step), the system wanders 
between trajectories, thus performing some sort of ran- 
dom walk in H with "repelling boundary conditions" at 
H m i n and " absorbing boundary conditions" on the walls 
(as negative densities are meaningless, the "death" of the 
system is declared when the trajectory hits the zero pop- 
ulation state for one of the species). This result was em- 
phasized by Gillespie ^l| for the important case where 
intrinsic stochastic fluctuations are induced by the dis- 
crete character of the reactants. In that case, the noise is 
multiplicative (proportional to the number of particles), 
and the system flows away from the center and eventu- 
ally hits one of the absorbing states at 0, or 0, 00. The 
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FIG. 3: The Lotka-Volterra phase space (left panel) admits a 
marginally stable fixed point surrounded by close trajectories 
(three of these are plotted). Each trajectory corresponds to 
single H defined in Eq. 1A2L where H increases monotoni- 
cally along the (dashed) line connecting the center with the 
a — wall, as shown in the lower right panel. In the up- 
per right panel, the period of a cycle T is plotted against H, 
and is shown to increase almost linearly from its initial value 
T — 2n/ y/JIa close to the center. 
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FIG. 4: The survival probability Q(t) is plotted versus time 
for a single patch noisy LV system. Eqs. lAH (with the sym- 
metric parameters) were integrated numerically (Euler inte- 
gration with time step 0.001), where the initial conditions are 
at the fixed point a = b = 1. At each time step, a small 
random number r/(t)At was added to each population den- 
sity, where r/(t) G [—A, A]. A typical phase space trajectory, 
for A = 0.5, is shown in the inset. The system "dies" when 
the trajectory hits the walls a = or b — 0. Using 300 dif- 
ferent noise histories, the survival probability is shown here 
for A = 0.5 (full line), A = 0.3 (dotted line) and A = 0.25 
(dashed line). Clearly, the survival probability decays expo- 
nentially at long times, Q(t) ~ exp(—t/r). As expected for a 
random walk with absorbing boundary conditions, 1 jr scales 
with A 2 . 



corresponding situation for a single patch Lotka-Volterra 
system with additive noise is demonstrated in Figure 4, 
where the survival probability Q(t) (the probability that 
a trajectory does not hit the absorbing walls within time 
t) is shown for different noise amplitudes. 



APPENDIX B: COUPLED OSCILLATORS 

The Lotka-Volterra system (Appendix A) is somewhat 
complicated, since the angular velocity depends not only 
on H, but also on the location along a trajectory. In 
order to simplify the discussion, let us introduce a toy 
model that imitates the main features of the real systems. 
Although that model does not allow for an absorbing 
state, we believe that it captures the basic mechanism for 
stabilization of spatially extended systems in the presence 
of noise. 

The toy model deals with the phase space behavior 
of diffusively coupled oscillators, where the angular fre- 
quency depends on the radius of oscillations. With addi- 
tive noise, the Langevin equations take the form: 



dx i 

~dt 

dx 2 

~5t 

dyi 
dt 

dy 2 
dt 



u(xi,yi)yi + Di(x 2 - Xt) + r)i(t) 
u(x 2 ,y 2 )y 2 + Dx(xi - x 2 ) + r/ 2 (t) 
-uj(x 1 ,y 1 )x 1 + D 2 (y 2 -y 1 )+T] 3 (t) (Bl) 
-u(x 2 ,y 2 )x 2 + D 2 {yi - y 2 ) + r)±(t). 



where all the 77-s are taken from the same distribu- 
tion. If the angular frequency is location independent, 
u)(x,y) — luq, the problem is reduced to coupled har- 
monic oscillators, a diagonalizable linear problem that 
admits two purely imaginary eigenvalues in the invari- 
ant, homogenous manifold. With the addition of noise, 
the random walk on that manifold is independent of the 
motion in the fast manifold, such that the radius of os- 
cillation diverges with the square root of time. As there 
are no "absorbing walls" here, the oscillation amplitude 
will grow indefinitely in the presence of noise. 

Now l et us de fine the oscillation radius for each patch, 
r i = V x i + Vi f° r * = 1)2) an d assume that the an- 
gular frequency depends only on that radius and is 9- 
independent [9i = arctg(yi/xi)]. With that, the total 
phase $ = 6\ + 9 2 decouples and the 3-dimensional phase 
space motion is dictated by the equations (we assume 
D X =D 2 = D and define d> = 61 - 9 2 ): 



r'i = -D(r x - r 2 cos(0)) + fji(t) 
r 2 = -D(r 2 - n cos(0)) + 772^) 



(B2) 



D 



nr 2 



siruj) — [co(r{) — w(r 2 )] + 



m _ m 

r\ r 2 



As before, all the 77-s are taken from the same distribu- 
tion. In the harmonic limit, w(ri) = uj{r 2 ) — ujq, the 
only noisy term for (f> vanishes at large r-s, and thus 
<f> approaches zero. Defining now the new coordinates, 
R = n + r 2 , r = n - r 2 , 



R= -2Dsin 2 (cj)/2)R + T] 
r = -2Dcos 2 ((f>/2)r + T], 



(B3) 
(B4) 
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FIG. 5: Histograms showing the probability to be at a dis- 
tance R from the origin as a function of R, for two coupled 
noisy oscillators, where u> = 1 + ar with D = 0.01 and vari- 
ous values of noise strength A and angular velocity gradient 
a. As expected, the phase space confinement is proportional 
to a, from a — 1 (triangles) to a = 0.5 (circles) to a = 0.1 
(solid line), all for the same level of noise A = 0.1. On the 
other hand, as predicted by the linear analysis close to the in- 
variant manifold, the confinement is noise independent, and 
the three solid lines corresponding to different levels of noise 
(A = 0.1, 0.5, 1 with the same a = 0.1 almost coincide. The 
inset shows the flow lines on the r = plane. The invariant 
manifold (f> — is stable, but there is a "maximally desynchro- 
nized" unstable orbit converging to the center at <f> — tt. If the 
expectation value of 2 deviates from zero, there is an effec- 
tive restoring force toward the center, and the noise-induced 
random walk on the = manifold is bounded. 

one notices that at the limit = the "restoring force" 
in the R direction vanishes and the phase space admits 
an attracting ld-manifold 4> = r = (the invariant man- 
ifold). The noise, thus, induces random walk on that R 
manifold, and since the walker cannot cross the origin, 
its displacement grows like \[t. 

In the generic case, however, where uj depends on r, 
the system behaves quite differently: if r\ ^ r-i, the 
two patches oscillate with different radial distance and 



desynchronize, i.e., (0 2 ) acquires, on average, some finite 
value. Migration, then, acts as a restoring force for an 
overdamped harmonic oscillator [Eq. IjBSJl ] and stabilizes 
the oscillations (if diffusion stabilizes the unstable fixed 
point this (noise independent) phenomenon is known as 
"oscillation death" |2fij). 

Following Eq. I|B3I) , close to the homogenous manifold 
(large diffusion, small noise limit), r 2 typical fluctuation 
around zero is of order A 2 /D. The oj{t\) — wlr?) term 
induces finite noise in the <j) equation of motion, and thus 
((f) 2 ) ~ (lu (r)) 2 A 2 /Z? 3 . This desynchronization, in turn, 
leads to the appearance of a finite restoring force on the 
invariant manifold R, as (4> 2 ) ^ 0, thus (far from the ori- 
gin) one finds (R 2 ) ~ D 2 /uj (r)) 2 . The small D instabil- 
ity (decoupled patches) manifest itself in the divergence 
of r 2 as D — > 0. It should be noted that, since both the 
restoring force and the noise in the invariant manifold 
are proportional to A 2 , the expected R distribution has 
to be noise independent at that limit, as demonstrated 
in Figures 5 and 2. All in all, if u> is radius dependent, 
the coupled oscillators' system stabilizes at some finite 
radius from the origin, giving rise to a soft "limit cycle" 
in the 4-dimensional phase space, as indicated in Figure 
5. The above considerations are valid only close to the 
homogenous manifold. For stronger noise, although the 
qualitative features of the system are the same, it was 
shown recently that |27| the effect of noise may shift the 
" intrinsic" frequency of the system, leading to some shift 
of the intrinsic frequency. 

Eqs. JHH may be generalized to include the case of an 
unstable focus (in order to imitate the dynamics of the 
Nicholson-Bailey |5( map) by adding a diagonal repulsive 
term to any variable (e.g., x\ — u>(xi,yi)yi + D%(x2 — 
x i) + Vi(t) + ea; i> where e measures the "strength" of 
the repulsion). The same analysis shows that, for small 
noise and small repulsion, a noise-induced transition will 
occur at e ~ (to 12 A 2 /D 2 ). If the noise is small enough, 
the desynchronization is weak and cannot stabilize a soft 
limit cycle, and thus the system is, for real populations, 
extinction-prone. Strong noise, conversely, stabilizes the 
system and ensures species conservation. 
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